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MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  THE 
DISCRETE  APPLICATION  OF  THE  AMSAA  GROWTH  MODEL 


1.  INTRODUCTION 

The  U.S.  Army  Materiel  Systems  Analysis  Activity  (AMSAA)  reliability 
growth  model  was  developed  under  the  leadership  of  Dr.  Larry  Crow  during  the 
mid  to  late  1970' s.  The  resultant  model  was  incorporated  into  Military 
Handbook  189,  Reliability  Growth  Management,  where  it  is  fully  explained. 

The  model  as  developed  by  Dr.  Crow  is  intended  for  use  with  continuous  data. 
However,  Dr.  Crow  later  developed  an  as  yet  unpublished  technique  for  applying 
the  model  to  discrete  data.  Briefly,  this  technique  assigns  a  sequence  number 
to  each  trial  of  the  test  and  uses  these  sequence  numbers  as  the  measures  of 
“time"  in  the  model.  The  growth  parameters  are  estimated  from  the  actual 
test  data  (the  number  of  test  intervals,  the  number  of  trials  in  each  interval, 
and  the  number  of  failures  in  each  interval)  by  the  maximum  likelihood  method. 

Mr  the  resultant  likelihood  equation  had  a  true  maximum,  then  any  of 
several  search  algorithms  could  be  applied  to  find  it.  However,  as  will  be 
shown  later,  the  likelihood  equation  can,  in  certain  regions,  be  unbounded. 

The  problem  only  becomes  solvable  if  the  maximum  over  a  restricted  region  is 
sought.  The  purpose  of  this  report  is  to  present  a  solution  to  the  problem 
and  a  computer  program  that  implements  that  solution. 

2.  THE  AMSAA  RELIABILITY  GROWTH  MODEL  AND  DISCRETE  DATA 

Use  of  the  AMSAA  model  assumes  the  test  is  divided  into  several  "intervals." 
Each  interval  represents  a  portion  of  the  test  conducted  on  fixed  configuration 
test  items.  One  interval  ends  and  another  begins  at  any  point  in  the  test 
where  significant  changes  are  made  in  the  test  items  in  an  attempt  to  change 
their  reliability  characteristics.  Therefore,  unless  changes  are  applied 
periodically,  it  is  inappropriate  to  use  the  model. 

The  model  itself  assumes  that  failures  occur  according  to  a  nonhomogeneous 
Poisson  process  with  a  Weibull  intensity  function  given  as: 

p(t)  =  XPtP’1  [1] 

where  A  is  the  scale  parameter  and  b  is  the  growth  or  shape  parameter.  From 
[1]  the  expected  number  of  failures  on  the  interval  [ Ti _ x  *  T i ]  is 
given  by: 

AT -  A  T j _ j B  [2] 

and  the  average  failure  rate  over  the  interval  is: 

A  ( Ti 6  -  TMB) 


Ti  -  Ti_i 


[3] 


To  apply  the  model  to  discrete  data,  the  sequence  numbers  of  the  last 
trial  in  each  interval  becomes  that  intervals  "end  time"  and  [3]  becomes  the 
average  probability  of  failure  for  each  trial  in  the  ith  interval.  If  K  is 
defined  as  the  number  of  intervals,  T-j  -  Tj_i  as  the  nunber  of  trials  in 
the  ith  interval  and  Mi  as  the  number  of  failures  in  the  ith  interval,  then 
a  likelihood  equation  for  the  test  can  be  constructed  using  [3]  as  the 
probability  of  failure  for  the  ith  interval  and  the  binomial  term  determined 
from  the  test  data  for  each  interval.  The  likelihood  function  then  is  given  by: 


The  problem  then  becomes  finding  the  values  of  x  and  p  which  maximize  [4]. 
To  simplify  the  mathematics  to  follow,  the  following  transforms  are  made: 


Define  X,  =  [5a] 

Define  N-j  =  T i - T -j _ i  [5b] 

Define  S,  =  t i - T i _ i - M j  [5c] 

Define  F-j  =  M-j  [5d] 


Then  [4]  becomes 


3.  MAXIMIZING  THE  LIKELIHOOD  FUNCTION 


Since  the  logarithm  of  F(X)  maximizes  at  the  same  point  as  F(X),  the 
values  of  x  and  p  which  maximize  [6]  will  also  maximize. 


K  K 

In*  =L=  v  (F-j)ln(xX-j)  +  V  (s^ln^  -  XX,) 
i=l  i=l 


In 


lnf^) Fi  -  1  n( ) Si 


[7] 


Note  that  the  last  term  is  a  constant  which  can  be  eliminated  without 
changing  the  values  of  X  and  p  which  will  maximize  [7].  Therefore,  the 
maximum  likelihood  estimates  of  X  and  p  (^,  $)  can  be  found  by  maximizing: 


K  K 

L  «  v  ( F j )  1  n ( X X i )  +  <•  (S-j)ln(N-j  -  XX,)  [8] 

i=l  i=l 
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.'ll 


**v 
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>1 
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To  find  x  and  8  from  [8]  would  require  that  [8]  be  bounded  and  continuous. 
However,  it  is  readily  seen  that  for  any  fixed  value  of  8.  if  X  increases 
without  bound,  the  point  is  eventually  reached  where  N^-xXi  =  0.  At  this 
point  L  becomes  discontinuous.  In  fact,  by  continuing  this  line  of  reasoning, 
it  may  be  shown  that  for  any  fixed  8  there  are  K  discontinuities  corresponding 

to  the  points  Nj  =  xXj,  N2  =  XX2 . N|<  =  XX|(.  To  be  able  to 

solve  [8]  for'ft  and  then  some  restriction  on  [8]  is  needed  that  is  not 
inconsistent  with  the  intent  of  the  solution.  As  it  turns  out,  a  real  world 
limitation  on  [8]  provides  exactly  such  a  restriction. 


From  [3],  [5a],  and  [5b],  it  is  noted  that  xX^/N^  is  an  estimate  of  the 
probability  of  failure  over  the  ith  interval.  Since  it  is  a  probability,  it 
can  be  restricted  to  (0,1)  without  loss  of  generality  for  our  problem  (indeed, 
doing  so  enhances  the  solution).  From  this  restriction,  limits  are  imposed 
on  x  namely. 


0  <  x  <  — 

Xi 


That  is,  x  must  be  such  that  it  is  less  than  N-j /X -j  for  each  i.  Coincidentally, 

this  is  precisely  the  restriction  required  to  eliminate  the  singularities  in 

[8].  So,  for  any  fixed  value  of  e.  there  is  an  upper  limit,  Xmax(8), 

imposed  on  x.  It  may  be  shown  that  xmax(e)  =  Nj/Xi  for  8<1 

and  Xmax(e)  *  Nk/X|<  for  B>1*  [Note  that  Nk/Xk  =  N]/Xi  =  1  =  xmax(l) 

at  8  =  1.]  This  yields  a  finite  function,  [8],  over  a  restricted  region 

defined  by  [9].  t  and  ^  can  now  be  sought. 


Recalling  that  the  X^'s  are  functions  of  8  and  defining  X  ^  to  be  the 
partial  of  Xj  with  respect  to  8,  the  first  partials  of  [8]  are  given  by: 


6L 

K 

y 

[ih 

K 

l 

i  =  l 

Si  xi 

6X 

i  =  l 

XXi 

N-j-xXj 

61 

K 

I 

i=l 

Fi*x‘i  _ 
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si  X  X* 
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XX  j 
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Continuing,  the  2nd  partials  of  [8]  are  given  by: 
*2i  f .  c. v  2 


62L  Fi 

-r  -  - 1  4  - 1 


(N.j -xX.j )‘ 


-  •  l 


FiTi6Ti-iB(lnTr1"Tn) 

Y  2 


xSiX  i 


x2si(x,i)2 


(N-xXi) 


(N-xXi )‘ 


=  .  V 


[14] 


6X68 


*S1*j  X'i  _  y  SiX'i 

(  N-XXj )  ^  (  N-XX-j ) 


The  two  terms  of  [13]  produced  by  taking  the  derivative  of  the  first  term  of 
[11]  have  been  combined  to  make  it  clear  that  the  sign  of  the  combined  term 
is  negative.  In  fact,  since  T  ^  p ,  Fj,  Sj,  X,  X  j.  X  j  and  (lnT^-lnT^ ) 
can  all  be  shown  to  always  be  positive  in  the  restricted  region,  it  may  be 
seen  that  [12],  [13],  and  [14]  are  all  negative  everywhere  in  the  restricted 
region.  This  indicates  that  [8]  is  concave  down  and  there  are  no  points  of 
inflection  over  the  entire  region.  Since  [8]  is  also  finite  and  continuous 
over  the  entire  region,  there  are  no  local  maxima  and,  if  there  is  a  point 
such  that  [10]  and  [11]  are  both  zero,  it  is  the  global  maximum.  These 
results  usually  indicate  the  use  of  nonlinear  programming  algorithms  to  find 
the  maximum.  However,  each  of  several  techniques  tried  had  problems  involving 
either  inordinate  computational  requirements  or  difficulty  in  dealing  with 
the  region  boundry.  Following  these  problems,  the  possibility  of  using  a 
bracketing  technique  was  investigated. 


Consider  the  nature  of  [10]  for  some  fixed  b»  say  e0-  Since  6  is 
held  constant,  the  X-j's  are  also  constant. 

If  x  is  set  to  zero,  [10]  becomes  +»  because  of  the  first  term.  If  X 
is  set  to  xmax(?0)then  either  N^XX^  =  0  or  N^-xXk  =  0  and  [10] 
becomes  Now  [12]  shows  that  [10]  monotonical ly  decreases  as  x  increases 

and  6  is  held  constant.  Therefore,  the  existence  of  some  x0  such  that 
[10]  evaluated  at  (x0,  B0)  equals  zero  is  guaranteed.  Since  this 
argument  can  be  made  for  any  b0>  it  indicates  that  there  exists  some  function 
of  6.  call  it  X*(b),  such  that  [10]  evaluated  at  ( X* ( p0 ) »  Bo) 
is  equal  to  zero.  Now,  if  x  >  x*(b0)  then  [10]  becomes  negative  (because 
it  decreases  monotonical ly)  and  L,  consequently,  decreases.  If  X  <  X*(Bo) 
then  [10]  becomes  positive  and  L  decreases.  That  is,  x*  (bo)  defines  the 
value  of  X  which  maximizes  [8]  for  8  =  B0-  The  global  maximum  of  [8]  over 
the  restricted  region  then  must  occur  at  one  of  these  points.  By  algebraic 
manipulation  of  [10]  and  [5a],  it  can  be  shown  that  Figure  1  depicts  the 
restricted  region  as  viewed  looking  down  on  the  x,  p  plane. 


4 


figure.  1. 


Figure  1  represents  the  general  shape  of  the  region.  Specific  values  shown 
must  occur.  That  is,  Xmax(0)=Ti, 


;**(0)  =  .  ,  xmax(l)  =  1  and 

^i 

X*(l)  =  .  always.  The  only  possible  deviation  from  Figure  1  for 

rNi 


the  specific  points  is  that  x*(0)  is  less  than  1  in  a  few  special  cases. 

For  most  legal  datasets  (that  is,  at  least  one  failure  and  at  least  one 
success)  A* ( 0)  is  between  1  and  T.  The  shape  of  Xmax(c)  can  be  verified 
by  simply  plotting  xmax(s)  =  N ^/ ( T ~ T K_ j e)  for  p<l  and  xmax(p) 

=  Ni/Tj6  for  8<1.  It  has  already  been  shown  that  X*(6)  must  exist 
in  the  open  interval  (0,  Xmax{e)).  Finally,  the  shape  of  X*(p)  can  be 
inferred  as  follows.  Select  and  So  and  X*(Bn).  At  this  point  [10] 
is  zero.  If  8  is  increased  to  80+<5  then  [10]  becomes  negative,  as 
indicated  by  [14]  being  negative  (that  is,  an  increase  in  8  causes  a  decrease 
in  [10]).  Now.  because  [12]  is  negative  everywhere,  x  has  to  be  decreased  to 
bring  [10]  back  to  zero.  Therefore,  X*(e+6)  <  X*(B)  for  positive  6. 


All  this  indicates  that  X*(P0)  can  be  found  by  a  binary  search 
technique.  It  is  known  that  [10]  is  +°°  at  x=0,  B  =  Bo  dnd 
at  x  =  Xmax(p0),  8=Po*  [10]  is  then  evaluated  at  x=xmax(p)/2 

and  it  is  positive,  then  X*(p0)  is  known  to  be  in  the  interval  (xnax( pc)/2 , 
Xir.ax(e0)).  Conversely,  if  [10]  is  negative  X*(p0)  is  in  the  interval 
(0,  Xmax(  ?jQ  )/2) .  This  can  be  continued  until  either  X*(e0)  is 
found  ([10]  equals  zero)  or  the  width  of  the  interval  known  to  contain  X*  ( p,0 ) 
is  less  than  the  accuracy  requirement  in  the  solution. 

This  bracketing  technique  provides  a  means  of  locating  x*(e)  for 
any  given  g.  What  remains  is  to  develop  a  legitimate  technique  for  searching 
along  the  function  X*(b)  for  the  point  that  maximizes  [8].  To  this 
end  ,  consider  Figure  2 . 


Figure  2. 


Assume  point  B  exists  as  a  point  on  x*(b)  such  that  [11]  is  negative. 

That  is,  at  B  the  slope  of  L  with  respect  to  X  is  zero  and  the  slope  with 
respect  to  p  is  negative.  Now  assume  that  point  A  exists  such  that  L  evaluated 
at  A  (I_a)  is  greater  than  Lb.  Since  there  are  no  points  of  inflection 
(and,  therefore,  no  saddle  points)  in  the  restricted  region,  there  must  exist 
some  point  C  between  A  and  B,  such  that  Lq  >  Lg»  at  any  arbitrarily  selected 
distance  from  B.  Figure  1  depicts  C  as  a  point  on  X*(b)  but,  for  this 
argument,  that  is  not  required.  Since  C  must  exist  at  any  radius  arbitrarily 
close  to  B,  it  must  exist  at  a  radius  such  that  point  B'  is  arbitrarily 
close  to  B.  Since  [13]  is  continuous  on  the  region  a  B'  can  always  be  selected 
so  that  [11]  is  still  negative.  Now  Lg'  <  Lb  by  definition  of  X*(b). 

Also,  since  [14]  is  negative,  increasing  8  from  6b  To  6c  causes  L  to 
decrease  (since  the  slope  of  L  with  respect  to  g  is  negative).  Therefore, 

Lg  <  Lg’  <  Lg.  Since  this  violates  the  requirement  that  Lq  >  Lg  the 
initial  assumption,  that  a  point  can  exist  at  ba  >  6b  when  [H]  is 
negative  and  La  >  Lg  leads  to  a  contradiction.  Using  a  similar  argument, 
it  can  be  shown  that  if  [11]  is  positive  at  some  point  (X*(Bg),  6b) 
then  there  cannot  exist  any  point  (x*(6a)»  6a)  where  Bg  >  6a 
and  La  >  Lb- 


As  an  end  result,  if  [11]  evaluated  at  some  (x*(p0),  e0)  is 
negative,  then  the  global  maximum  will  occur  at  some  p  <  ?0.  Conversely, 
if  [11]  evaluated  at  (x*(60) .  So)  is  positive  then  any  global 
maximum  that  exists  will  occur  at  some  B>S0-  Note  that  there  is  no 
guarantee  that  a  point  exists  such  that  [10]  and  [11]  are  both  zero.  It  has 
not  been  proved  that  [11]  cannot  be  asymptotic  to  zero  as  e  increases  without 
bound.  In  this  case,  however,  there  is  a  realistic  limit  imposed  on  b  during 
the  search.  This  limit  is  discussed  in  the  next  section. 

This  completes  the  proof  that,  theoretically,  X,  s'  can  be  found  by  a 
two  dimensional  binary  search  technique.  The  next  section  of  this  paper 
deals  with  implementing  this  theoretical  approach  on  a  digital  computer. 

4.  FINDING  X  AND  ?  ON  A  DIGITAL  COMPUTER 

The  previous  section  showed  that,  in  theory,  [8]  has  a  maximum  which 
can  be  located  using  a  two  dimensional  search  technique.  This  theory  in  its 
pure  form,  however,  depends  on  being  able  to  locate  for  any  e  a  corresponding 
X  such  that  [10]  is  negative,  but  finite,  and  another  such  that  [10]  is 
positive,  but  finite.  It  also  depends  on  being  able  to  find  a  value  of  £  such 
that  [11]  evaluated  at  (x*(b),  B)  is  negative.  The  existence  of 
the  first  two  points  was  proved  in  the  previous  section.  It  is  possible, 
however,  that  while  the  points  exist  it  may  not  be  possible  to  compute  them 
on  a  digital  computer.  It  is  possible  to  construct  a  problem  such  that  for 
some  e  x*(e)  is  so  close  to  Xmax(g)  as  to  be  beyona  the  accuracy 
of  the  computer  to  discriminate  between  them.  The  theory,  for  instance,  is 
satisfied  by  Xmax(e)  =  XQ  and  X* ( 6 )  =  XQ  -  2"  .  (Assume  X0 

to  be  on  the  order  of  1.)  On  a  digital  computer  with,  say,  23  bits  for  the 
mantissa  of  a  floating  point  variable,  the  actual  value  computed  for  xmax(f) 
may  be  as  small  \Q  -  2"  ,  which  is  less  than  x0  -  2“  .  That  is,  because 

of  truncation  error  involved  in  the  calculation  of  Xmax(e),  it  is 
possible  that  [10]  evaluated  at  (xmax(p),  $'  would  be  positive 
rather  than  negative.  This  result  indicates  that  x*(e)  is  too  close 
to  xmax(p)  for  the  computer  to  discriminate  between  the  two. 

Similarly,  it  may  be  that  x*(e)  is  so  close  to  zero  so  as  to  be 
indistinguishable  from  zero  by  a  computer.  Based  on  experience  with  actual 
test  cases,  it  is  most  likely  that  X*(b)  being  computationally  equal 
to  zero  will  be  encountered  when  the  actual  solution  for  b  is  very  large 
(say  e  >  3).  This  will  cause  Xmax(B),  and  consequently  x*(p), 
to  pass  very  close  to  the  B  axis. 

Finally,  the  existance  of  a  point  on  x*(b)  such  that  [11]  is 
negative  has  not  been  proved. 

Fortunately,  all  of  these  situations  have  realistic  work  arounds.  In 
the  first  case,  the  solution  is  to  set  X*(b)  =  Xmax(g)  as  computed 
and  continue  the  algorithm.  Because  [8]  is  so  very  well  behaved  in  the 
restricted  region,  small  errors  in  the  location  of  any  particular  parameter 
estimate  has  a  correspondingly  small  effect  on  the  solution.  Experience  with 
actual  data  indicates  that  truncation  error  in  calculating  Xmax(B)  is 
enough  to  prevent  [10]  from  evaluating  to  In  the  second  and  third  cases, 

an  unsolvable  problem  is  probably  trivial.  In  several  years  of  experience 
with  the  program  outlined  later  there  has  never  been  a  case  where  x*(f) 


became  computationally  zero.  Most  of  these  runs  have  been  made  with  a  cutoff 
of  2  for  8  •  That  is,  if  [11]  was  positive  at  (x*(2),  2)  the  program  was 
stopped  and  a  message  that  ^  was  greater  than  2  was  printed.  In  other  words, 
in  every  case  ever  run  x*(2)  was  non-zero.  In  real  world  terms,  if  6  >  2 
then  the  actual  value  really  does  not  matter.  If  growth  is  occurring,  then  p  < 
For^  >  1  the  reliability  of  the  item  is  actually  degrading  as  "improvements" 
are  made.  In  short,  determining  that  £  >  2  should  be  a  sufficient  result 
without  knowing  the  actual  value  of  $  . 

5.  A  COMPUTER  PROGRAM  FOR  COMPUTING  ^  AND  £ 

The  appendix  contains  a  computer  program  to  compute  x  and  8-  The 
program  is  written  in  FORTRAN  77  and  was  developed  on  a  VAX  11/780  under  the 
VMS  operating  system. 

Main  Program  -  The  main  program  is  provided  for  convenience  and  provides 
no  part  of  the  solution.  It  may  be  freely  replaced  with  any  other  program 
which  sets  up  the  vectors  N  and  F,  provides  for  vectors  X  and  DX  and  then 
calls  SUBROUTINE  MLHE. 

SUBROUTINE  MLHE  -  This  is  the  routine  that  drives  all  others  involved  in 
the  solution.  Statements  30  to  35  initialize  variables  LAMBDA,  BETA,  EPS, 

TF,  and  TN. 

Statements  36  to  47  do  some  error  checking  on  the  input  vectors,  N  and  F. 
An  error  condition  (MSG*0)  is  returned  if  1)  any  interval  has  less  than  one 
trial,  2)  any  interval  has  less  than  zero  failures,  3)  any  interval  has  more 
failures  than  trials,  4)  all  trials  were  failures,  or  5)  no  trials  were 
failures.  Any  of  these  conditions  will  cause  the  computation  of  x  and  8  to 
not  take  place. 

Statements  48  -  53  set  BETA  to  2  (the  artificial  upper  limit  discussed 
in  the  last  section)  and  computes  the  corresponding  x*  as  well  as  the  value 

of  [11]  by  calling  SUBROUTINE  GETLAM.  If  the  value  of  [11]  (variable  DB)  is 

positive  it  indicates  that  ^  >  2.  In  this  case,  the  routine  is  terminated 
with  MSG  =  3,  indicating  that  £  is  greater  than  2. 

Statements  54  -  69  perform  the  binary  search  for  Beginning  with 
BH=2  and  BL=zero.  The  value  of  DB  is  determined  at  (BH+BL)/2.  If  DB  is 

greater  than  zero,  BL  is  set  to  the  average.  If  DB  is  less  than  zero,  BH  is 

set  to  the  average.  The  loop  is  then  executed  again  using  the  new  BH  and  BL. 
The  loop  is  terminated  on  either  of  two  conditions:  1)  The  difference 
between  BH  and  BL  is  less  than  the  accuracy  requirement  (EPS),  or  2)  The 
number  of  trips  through  the  loop  exceeds  100.  The  second  condition  is 
necessary  to  prevent  infinite  loops  under  certain  conditions.  If  too  much 
accuracy  (variable  NSIG)  is  requested,  it  is  possible  to  get  into  the  situation 
where  BH-BL  computes  to  a  value  greater  than  EPS,  but  (BH+BL)/2  computes 
to  BH  or  BL.  This  situation  occurred  on  a  VAX  11/780  under  VMS  FORTRAN  when 
NSIG  was  set  to  6. 


Statements  70  -  71  take  the  average  of  the  last  values  of  BH  and  BL  and 
call  GETLAM  to  compute  LAMBDA.  This  represents  the  solution  to  the  problem 
and  a  return  to  the  main  program  is  taken. 

SUBROUTINE  GETLAM  -  Statements  5-10  are  executed  if  BETA  *  0.  This 
can  occur  when  NSIG  is  too  large  for  the  machine  and  EPS  computes  to  zero. 

In  this  case,  it  is  easily  shown  that  X*  is  given  by  the  expression  in 
statement  9.  At  8=0,  all  the  X-j ' s  become  zero  except  for  Xi,  which  becomes  1. 
[10]  is  easily  solved  for  x.  Similarly,  since  all  Xf  =  0  it  is  obvious 
that  [11]  becomes  +•>.  Since  the  search  algorithm  only  depends  on  the  sign 
of  DB  and  DL  and  not  their  magnitude,  DB  is  set  to  an  artificial  positive 
value  and  a  return  is  taken  to  MLHE. 

,  Statements  11  -  22  compute  the  values  of  the  X^'s  (vecter  X)  and  the 
X  jS  (vector  DX)  for  the  given  value  of  BETA. 

Statements  23  -  28  handle  the  situation  where  8=1.  In  this  case,  X^  *  N, 
for  all  i  and  x*  =  total  failures/total  trials  Is  easily  shown  by  setting 
[10]  to  zero,  setting  X^  =  N^  and  solving  for  x. 

Statements  29  -  51  perform  the  same  binary  search  algorithm  on  x  that 
is  performed  on  8  by  MLHE.  The  loop  termination  conditions  are  the  same  and 
for  the  same  reasons.  Since  It  can  be  shown,  mathematically,  that  [10]  =  +« 
at  x=0  and  [10]  =  —  at  x*Xmax,  these  values  are  selected  as  the 
initial  limits  on  the  bracket.  Statement  54  calls  value  to  calculate  DB 
using  the  given  BETA  and  LAMBDA  *  x*. 

SUBROUTINE  VALUE  -  The  purpose  of  this  routine  Is  to  compute  the  values 
of  the  first  partials  ([10]  and  [11])  at  LAMBDA  and  BETA.  Actually,  a  scaled 
value  of  DB  is  calculated.  This  is  done  to  exploit  commonality  between  [10] 
and  [11].  Since  the  search  algorithms  depend  only  on  the  sign  of  the 
derivatives,  this  scaling  does  not  effect  the  solution.  To  understand  the 
scaling,  consider  [10]  and  [11]  rewritten  as  follows: 


Note  that  the  terms  inside  the  parenthesis  of  both  [15]  and  [16]  are  Identical. 
Both  equations  can  be  computed  with  the  same  algorithm  by  just  substituting 
Xi  or  xX  <  In  the  first  position  of  each  term.  Note  also  that,  since  x  >  0, 
dividing  [16]  by  x  will  not  effect  whether  It  evaluates  to  a  negative,  zero 
or  positive  value.  This  1$  exploited  by  VALUE  which  uses  a  dummy  vector  (C) 
to  toggle  between  X^  and  X  ^  depending  on  whether  DL  or  DB  Is  desired.  To 
compute  DL,  VALUE  Is  called  with  C  ■  X.  To  compute  the  scaled  value  of  DB 
VALUE  is  called  with  C  «  DX. 
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SUBROUTINE  MESSAGE  -  This  subroutine  adds  nothing  to  the  solution  of  X 
and  b  .  It  is  provided  for  convenience  and  to  explain  the  meanings  of  the 
various  values  of  MSG.  It  may  be  freely  replaced.  It  is  invoked  only  from 
the  provided  optional  main  program. 


10 


’  *'i 


UE 

,<N 


APPENDIX 


>  .t 

rv:? 

r> 


if 


The  next  page  is  blank. 


11 


PROGRAM  MLHEPRO 
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C 

C  THIS  PROGRAM  COMPUTES  THE  MAXIMUM  LIKLIHOOD  ESTIMATES  OF  THE 
C  PARAMETERS  FOR  THE  AMSAA  RELIABILITY  GROWTH  MODEL  AS  APPLIED  TO 
C  DISCRETE  DATA.  THE  VARIABLES  USED  ARE  DEFINED  AS  FOLLOWS: 

C 

C  IOPT— INPUT— AN  INTEGER  VALUE  WHICH  SELECTS  ONE  OF  TWO  OPTIONS. 

C  IOPT=0  SELECTS  THE  SPECIAL  CASE  FOR  WHICH  ALL 

C  INTERVALS  CONTAIN  EXACTLY  ONE  TRIAL.  IOPT«l 

C  SELECTS  THE  GENERAL  CASE  FOR  WHICH  EACH  INTERVAL 

C  CAN  CONTAIN  ANY  NUMBER  OF  TRIALS. 

C  K-INPUT-AN  INTEGER  VALUE  WHICH  SPECIFIES  THE  NUMBER  OF 
C  INTERVALS  IN  THE  CURRENT  DATA  SET. 

C  N-INPUT-A  REAL  ARRAY  WHICH  CONTAINS  THE  NUMBER  OF  TRIALS 
C  IN  EACH  INTERVAL.  THAT  IS,  N(I)  IS  THE  NUMBER 

C  OF  TRIALS  IN  THE  ITH  INTERVAL. 

C  F-INPUT-A  REAL  ARRAY  WHICH  CONTAINS  THE  NUMBER  OF  FAILURES 
C  IN  EACH  INTERVAL.  THAT  IS,  F(I)  IS  THE  NUMBER  OF 

C  TRIALS  IN  THE  ITH  INTERVAL. 

C 

C  NOTE  THAT  THE  MAXIMUM  NUMBER  OF  INTERVALS  ALLOWED  IS  500.  TO 
C  INCREASE  THE  NUMBER  ALLOWABLE  JUST  RESET  THE  DIMENSION  LIMITS 
C  ON  EACH  OF  THE  ARRAYS  IN  THE  DIMENSION  STATEMENT. 

C 

C  THE  AMSAA  RELIABILITY  GROWTH  MODEL  FOR  DESCRETE  DATA  WAS  DEVELOPED 
C  BY  DR.  LARRY  CROW  OF  AMSAA. 

C 

C  THE  COMPUTER  PROGRAM  TO  SOLVE  THE  MODEL  WAS  DEVELOPED  BY  MR.  WILLIAM 
C  CLAY  OF  AMSAA.  ANY  QUESTIONS  CONCERNING  THE  PROGRAM  SHOULD  BE 
C  DIRECTED  TO  MR.  CLAY  AT  AV298-6887  OR  COMMERCIAL  (301)  278-6887 

C 

C 

DIMENSION  N ( 500 ) ,F<500) ,S<500> ,X(500> ,DX(500> 

CHARACTER* 1  ANS 
REAL  N,  LAMBDA 
NSET  =  0 
2  PRINT  1 

1  FORMAT (///'  ENTER  1  FOR  GENERAL  CASE  MODEL', 

+  / '  ENTER  0  FOR  ROUND-BY-ROUND  MODEL ' , 

+  /'  ENTER  ANYTHING  ELSE  TO  STOP',/) 

C 

C  READ  IN  THE  OPTION  AND  THE  NUMBER  OF  INTERVALS  FOR  THE  NEXT 
C  DATA  SET. 

C 

READ  * , IOPT 

IF  (IOPT  .LT.  0  .OR.  IOPT  .GT.  i)  STOP 
PRINT  110 
C 

C  NOTE:  THE  *  IN  THE  FORMAT  STATEMENT  CAUSES  THE  NEXT  READ  TO  BE 
C  SOLICITED  ON  THE  SAME  LINE  AS  THE  PRINT.  IF  THIS  OPTION 

C  IS  NOT  IMPLEMENTED  ON  YOUR  COMPUTER,  THE  *'S  SHOULD  BE 

C  REMOVED  FROM  THE  FORMAT  STATEMENTS. 

C 

110  FORMAT (/'  ENTER  THE  NUMBER  OF  INTERVALS  -  '*) 
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055  READ  *,  K 

056  IF  <K  .LT.  2)  THEN 

057  PRINT  120 

058  120  FORMAT ('  MODEL  REQUIRES  AT  LEAST  2  INTERVALS') 

059  GO  TO  2 

060  END IF 

061  NSET  =  NSET  +  1 

062  IF  ( I OPT  .EQ.  1)  THEN 

063  C 


064 

C 

IF 

GENERAL  CASE  OPTION  MAS  SELECTED  (IOPT  *  1).  READ  IN  THE 

065 

C 

NUMBER  OF  TRIALS  AND  THE  NUMBER  OF  FAILURES  FOR  EACH  INTERVAL 

066 

C 

IN 

ORDER. 

067 

c 

068 

DO  6  I  =  1,K 

069 

PRINT  15,  I 

070 

15 

FORMAT ( '  ENTER  NO.  ROUNDS,  FAILURES  FOR  INTERVAL  ',13,' 

071 

READ  *,  N ( I ) ,  F ( I ) 

072 

S(I)  =  N ( I )  -  F ( I ) 

073 

6 

CONTINUE 

074 

PRINT  3, NSET, K 

075 

3 

FORMAT!/'  DATA  SET  ',15,'  NUMBER  OF  GROUPS  IS  ',15,/, 

076 

« 

1H  , 3 ( '  N  S  F ' , 7X ) ) 

077 

PRINT  4, (N(I) ,S(I) ,F(I) ,I=1,K) 

078 

4 

FORMAT ( 1H  , 3F6. 0 , 7X , 3F6. 0 , 7X , 3F6. 0> 

079 

ELSE 

030 

c 

081 

c 

FOR 

THE  SPECIAL  CASE  OPTION  (IOPT  =  0)  READ  IN  THE  NUMBER  OF 

082 

c 

FAILURES  FOR  EACH  INTERVAL-  IT  IS  ASSUMED  THAT  THE  NUMBER  OF 

083 

c 

TRIALS  FOR  EACH  INTERVAL  IS  EXACTLY  ONE. 

084 

c 

085 

DO  199  I  =  1,  K 

086 

F  ( I )  =  0. 

087 

199 

CONTINUE 

088 

200 

PRINT  201 

089 

201 

FORMAT ('  ENTER  SEQUENCE  NUMBER  OF  FAILED  ROUND  (0  TO  QUIT) 

090 

READ  * ,  INDF 

091 

IF  (INDF  .LE.  0)  THEN 

092 

GO  TO  203 

093 

ELSE  IF  (INDF  .LE.  K)  THEN 

094 

F(INDF)  =  1. 

095 

ELSE 

096 

PRINT  202 

097 

202 

FORMAT ( '  INDEX  NUMBER  EXCEEDS  NUMBER  OF  ROUNDS') 

098 

END  IF 

099 

GO  TO  200 

100 

203 

CONTINUE 

101 

DO  11  I  =  1  ,K 

102 

N  (  I >  =  1. 

103 

11 

CONTINUE 

104 

PRINT  12, NSET, K 

105 

12 

FORMAT (/'  DATA  SET  ',15,'  NUMBER  OF  ROUNDS  IS  ',15) 

106 

PRINT  13, <F(I> ,1=1, K) 

107 

13 

FORMAT (20F4.0) 

108 

END  IF 

109 

CALL  MLHE (LAMBDA, BETA, K,N,F,X,DX, 4, DB,DL, MSG) 

14 


IF  (MSG  .EQ.  0)  THEN 
PRINT  5, LAMBDA, BETA 

FORMAT (/ '  ESTIMATED  LAMBDA  =  *,F5.3,'  ESTIMATED 
PLAST  =  1 .  -  LAMBDA  *  X  <K)  /  NOO 
PRINT  7, PL AST 

FORMAT ('  FINAL  ESTIMATE  OF  PROBABILITY  OF  SUCCESS 
ELSE 

CALL  MESSAG (MSG) 

END IF  i 

GO  TO  2 


ESTIMATED  BETA  «  ' ,F5. 3) 


' ,F10.8) 


SUBROUT I NE  MLHE (LAMBDA , BETA , K ,N,F, X ,DX ,NSIG,DB,DL, MSG ) 
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C 

C  THIS  ROUTINE  COMPUTES  THE  MAXIMUM  LIKLIHOOD  ESTIMATES  OF 
C  THE  TWO  PARAMETERS,  LAMBDA  AND  BETA.  VARIABLES  ARE  AS  FOLLOWS: 

C 

C  LAMBDA-OUTPUT-THE  MAXIMUM  LIKILIHOOD  ESTIMATE  OF  LAMBDA. 

C  BETA-OUTPUT-THE  MAXIMUM  LIKLIHOOD  ESTIMATE  OF  BETA. 

C  K-INPUT-THE  NUMBER  OF  INTERVALS. 

C  N-INPUT-AN  ARRAY  CONTAINING  THE  NUMBER  OF  TRIALS  FOR  EACH  INTERVAL 

C  F-INPUT-AN  ARRAY  CONTAINING  THE  NUMBER  OF  FAILURES  FOR 

C  EACH  INTERVAL. 

C  X,DX -WORKING-ARRAYS  OF  AT  LEAST  K  WORDS  EACH  USED  AS  WORKING 

C  STORAGE  WITHIN  THE  SUBROUTINE. 

C  NSIG— INPUT-THE  NUMBER  OF  SIGNIFICANT  DIGITS  REQUIRED  IN 
C  THE  SOLUTION.  NOTE:  COMPUTATIONAL  TIME  AND  TRUNCATION 

C  PROBLEMS  INCREASE  SIGNIFICANTLY  IF  NSIG  IS  MUCH  GREATER 

C  THAN  4. 

C  DB-OUTPUT-A  SCALED  VALUE  OF  THE  PARTIAL  OF  L  WITH  RESPECT  TO  BETA 

C  THE  VALUE  IS  NOT  CORRECT  (UNLESS  IT  IS  ZERO) ,  BUT  THE 

C  SIGN  IS. 

C  DL-OUTPUT-THE  PARTIAL  OF  L  WITH  RESPECT  TO  LAMBDA.  VALUE  AND  SIGN 

C  ARE  CORRECT. 

C  MSG-OUTPUT-AN  INTEGER  VALUE  DESCRIBING  THE  CONDITION  OF 

C  TERMINATION.  SEE  SUBROUTINE  MESSAG  FOR  MORE 

C  DETAIL. 

C 

REAL  LAMBDA ,N 
LOGICAL  LAST 

DIMENSION  N(K) ,F(K) ,X (K) ,DX(K> 

LAMBDA  =  0. 

BETA  =  0. 

MSG  =  0 

EPS  =  10.  **  (-NSIG)  /  2. 

TF  =  0. 

TN  =  0. 

DO  1  I  =  1 , K 

IF  (N ( I >  .LE.  0. )  MSG  =  1 
IF  (F ( I )  .LT.  0. )  MSG  =  5 
IF  (F ( I )  .GT.  N ( I ) )  MSG  =  4 
IF  (MSG  .NE.  0)  RETURN 
TF  =  TF  +  F ( I ) 

TN  =  TN  +  N ( I ) 

1  CONTINUE 

IF  (TF  .EQ.  0.  .OR.  TN  .EQ.  TF)  THEN 
MSG  =  2 
RETURN 
END  IF 
BETA  =  2. 

CALL  GETLAM (LAMBDA, BETA, K,N,F,X,DX, EPS, DB,DL,TF,TN) 

IF  (DB  .GE.  0.)  THEN 
MSG  =  3 
RETURN 
END  IF 
BH  *  2. 


TRUE 


LAST  =  .FALSE. 

LOOPS  =  0 

2  if  (BH-BL  .LE.  EPS  .OR.  LOOPS  .GE.  100)  LAST  *  . 
LOOPS  =  LOOPS  +  1 
BETA  =  <BH  +  BL)  /  2. 

CALL  GETLAM ( LAMBDA , BET A , K , N , F , X , DX , EPS , DB , DL , TF  t  TN ) 
IF  <DB. EQ. 0. >  THEN 
RETURN 

ELSE IF  (DB  .GT.  0.)  THEN 
BL  =  BETA 
ELSE 

BH  =  BETA 
END  IF 

IF  (.NOT.  LAST)  GO  TO  2 
10  BETA  *  (BH  +  BL)  /  2. 

CALL  GETLAM ( LAMBDA , BETA ,  K ,  N ,  F  ,  X  ,  DX  , EPS , DB f DL , TF , TN ) 
RETURN 
El 


001 

SUBRCHJT I NE  GETLAM  ( LAMBDA , BETA ,K,N,F, 

002 

REAL  LAMBDA ,N, LMIN, LMAX 

003 

LOGICAL  LAST 

004 

DIMENSION  N<K> ,F(K) ,X <K> ,DX (K> 

005 

IF  (BETA  .EQ.  0.)  THEN 

006 

LAMBDA  =  TF  *  N ( 1 )  /  ( TF  +  N ( 1 )  -  1 

007 

DL  “  0. 

008 

DB  «  1. 

009 

RETURN 

010 

END  IF 

011 

SUM  »  0. 

012 

PEX  *  0. 

013 

PDEX  *  0. 

014 

DO  1  I  *  1  ,K 

015 

SUM  =  SUM  ♦  N(I) 

016 

EX  *  SUM  **  BETA 

017 

DEX  =  EX  *  ALOG ( SUM ) 

018 

X(I)  =  EX  -  PEX 

019 

DXU)  =  DEX  -  PDEX 

020 

PEX  =  EX 

021 

PDEX  =  DEX 

022 

1 

CONTINUE 

023 

IF  (BETA  .EQ.  1.)  THEN 

024 

LAMBDA  =  TF  /  TN 

025 

DL  =  0. 

026 

CALL  VALUE  < LAMBDA, K,N,F, X ,DX ,DB) 

027 

RETURN 

028 

END  IF 

028 

LMIN  =  0. 

030 

LMAX  =  N ( 1 )  /  X(l) 

031 

IF  (BETA  .GT.  1.)  LMAX  *  N(K)  /  X(K> 

032 

UPRLIM  =  LMAX 

033 

LAST  «  .FALSE. 

034 

LOOPS  =  0 

035 

2 

IF  (LMAX  -  LMIN  .LE.  EPS  .OR.  LOOPS 

036 

LOOPS  =  LOOPS  +  1 

037 

LAMBDA  =  (LMAX  +  LMIN)  /  2. 

038 

CALL  VALUE ( LAMBDA, K , N, F, X , X ,DL> 

039 

IF (DL. EQ. 0. )  THEN 

040 

GO  TO  100 

041 

ELSE IF  (DL  .GT.  0.)  THEN 

042 

LMIN  «  LAMBDA 

043 

ELSE 

044 

LMAX  =  LAMBDA 

045 

END  IF 

046 

IF (.NOT.  LAST)  GO  TO  2 

047 

IF  (LMAX  .NE.  UPRLIM  .AND.  LMIN  .NE. 

048 

LAMBDA  *  (LMIN  +  LMAX)  /  2. 

049 

CALL  VALUE (LAMBDA, K,N,F,X,X,DL) 

050 

END  IF 

051 

100 

CALL  VALUE ( LAMBDA ,K,N,F,X,DX,DB) 

052 

RETURN 

053 

END 

.BE.  100)  LAST  =  .TRUE. 


0.)  THEN 


IB 


001  SUBROUTINE  VALUE (LAMBDA, K ,N,F, X ,C, DIRIV) 

002  REAL  LAMBDA, N 

003  DIMENSION  N(K) ,F(K> ,X(K) ,C(K) 

004  DIRIV  =  0. 

005  DO  100  I  =  1,K 

006  EX  =  LAMBDA  *  X  < I ) 

007  DIRIV  =  DIRIV  +  C(I>  *  N ( I )  *  (F(I)  -  EX)  /  ((NCI)  -  EX)  *  EX) 

008  100  CONTINUE 

009  RETURN 

010  END 
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SUBROUTINE  MESSAG(M) 

IF  <M. EQ. 1 )  PRINT  1 

1  FORMAT <‘  AN  INTERVAL  WAS  DETECTED  WITH  LESS  THAN  1  TRIAL' 
+  /'  NO  SOLUTION  WAS  ATTEMPTED') 

IF (M. EQ. 2)  PRINT  2 

2  FORMAT <'  TOTAL  NUMBER  OF  FAILURES  IS  ZERO  OR  TOTAL  NUMBER 

OF  SUCCESSES  IS  ZERO.  NO  SOLUTION  WAS  ATTEMPTED') 

IF (M. EQ. 3)  PRINT  3 

3  FORMAT ( '  THE  ESTIMATE  OF  BETA  EXCEEDS  2', 

+  /'  NO  SOLUTION  WAS  ATTEMPTED*) 

IF  <M. EQ. 4)  PRINT  4 

4  FORMAT ('  AN  INTERVAL  WAS  DETECTED  WITH  MORE  FAILURES  THAN 
♦TRIALS'/'  NO  SOLUTION  WAS  ATTEMPTED ' ) 

IF (M. EQ. 5)  PRINT  5 

5  FORMAT ('  AN  INTERVAL  WAS  DETECTED  WITH  LESS  THAN  ZERO  *, 
♦'FAILURES'/'  NO  SOLUTION  WAS  ATTEMPTED*) 

RETURN 

END 
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